Synchronization of weakly perturbed Markov chain oscillators 
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Rate processes are simple and analytically tractable models for many dynamical systems that 
switch stochastically between a discrete set of quasi stationary states; however, they may also 
approximate continuous processes by coarse-grained, symbolic dynamics. In contrast to limit-cycle 
oscillators that are weakly perturbed by noise, in such systems, stochasticity may be strong, and 
topologies more complicated than a circle can be considered. Here, we apply a second-order time- 
dependent perturbation theory to derive expressions for the mean frequency and phase diffusion 
constant of discrete-state oscillators coupled or driven through weakly time-dependent transition 
rates. We also describe a method of global control to optimize the response of the mean frequency 
in complex transition networks. 
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I. INTRODUCTION 



The emergence of an oscillating mean field in large en- 
sembles of noisy or nonidentical oscillators is considered 
the hallmark of synchronization as a collective effect 
However, the term synchronization is used differ- 
ently in different contexts. It may also refer to anything 
from statistical correlations to phase synchronization or 
complete synchronization. For individual oscillators, it is 
appropriate to define synchronization as the adjustment 
of frequency due to an interaction Q. Here it is in 
this sense that we will study synchronization for weakly 
driven or coupled Markov chains. If a system with 
frequency u/q is driven at frequency U3\, two effects may 
be observed: correlation of the state of the oscillator 
with the phase of the driving signal and adaptation 
of mean frequency f2. For deterministic limit-cycle 
oscillators, these effects correspond to phase and fre- 
quency synchronization and for coupling stronger than a 
critical value, phase and frequency locking are observed. 
These extreme cases are not observed in the presence 
of noise, although the Kramers theory for weak noise 
predicts an exponentially small frequency difference 
within the synchronization regions (Sec JIV C[) . On the 
other hand, strong fluctuations, such as those present 
in noise-induced oscillations and stochastic cycles over a 
few states or with heterogeneous transition rates, should 
be considered an integral part of the system. In this 
case, the first-order perturbation theory can predict 
the correlations induced by weak external perturbation. 
In addition, with the second-order perturbation theory 
presented in this paper, it is possible to quantify the 
change of frequency. 

Stochastic processes over discrete and finite sets of 
states have long been used as conceptual models for 
stochastic oscillators and proven to capture the essential 
mechanisms of synchronization. Discrete-state Markov 
rate models arc typically applicable for molecular ma- 
chines. Enzymes and motor proteins can be considered 



as molecular machines that undergo operation cycles 
consuming energy and performing work. Often these 
operation cycles are well modeled by a finite set of 
configurations and reaction rates, which describe the 
speed of transitions between the states. The transi- 
tions occur randomly whenever a substrate molecule 
binds, ATP is converted or through thermal activation. 
Molecular machines operating under non-equilibrium 
conditions are therefore stochastic oscillators that can 
be approximated by a continuous-time Markov chain 
model without detailed balance, i.e., with an average 
flow along the operation cycle (see Fig. QJ,). The mean 
frequency directly equates to the productivity of a 
molecular machine or the turnover rate of an enzymatic 
reaction, which can be increased or decreased depending 
on the situation. This can be achieved by driving the 
system purposefully at a resonant frequency. 

If the frequencies of two oscillators are sufficiently 
similar, mutual coupling can increase the coherence of 
the stochastic oscillations. This nontrivial collective 
effect, where the order is an emergent property of 
the coupled system, is an important mechanism to 
reduce noise in biological oscillators on mcso- and 
micro-scales. It has, for instance, been observed in 
ensembles of beating heart cells Q. In [tJ it was 
shown numerically that by coupling a few noisy gene 
regulatory circadian oscillators, phase diffusion for each 
oscillator decreases significantly. Another recent paper 
Q concludes that in some circumstances, the circadian 
protein phosphorylation cycle in cyanobacteria can 
dramatically enhance the robustness of the coexisting 
gene regulatory transcription-translation cycle. Both 
cycles were described by their respective rate processes, 
i.e., coupled stochastic oscillators of different types. 

In this paper we apply the second-order perturba- 
tion theory to the master equation of a Markov rate 
process to derive expressions for its mean frequency 
£1. The adjustment of this frequency to frequency ui\ 
of a deterministic or stochastic driving signal is clearly 
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FIG. 1. (Color online) Schematics of (a) a stochastic oscillator 
with three metastable states in a tilted harmonic potential. 
The transitions are Kramers rate processes and can be per- 
turbed through a change in the activation energy or noise 
strength. In this example the forward transitions are per- 
turbed in form of a traveling wave, (b) Nonsequential stochas- 
tic oscillator with four states. The diagrams on the right show 
possible lifts of the Markov chain to periodic lattices and asso- 
ciated Poincare sections. The red Poincare section (lower di- 
agram) disregards sub-threshold oscillations 1 — > 2 — >• 4 — > 1. 



identified as a parametric resonance phenomenon. For a 
simple jump process on a ring with discrete rotational 
symmetry we can also derive expressions for the phase 
diffusion constant in the case of external driving or 
mutual coupling between two stochastic oscillators. 
By exploiting the special structure of the second-order 
perturbation terms, we also present an algebraic method 
to maximize or minimize these terms under linear con- 
straints on the perturbation and quadratic constraints 
on its power. Thus, it is possible to optimize pump 
flows in a diffusion process over a directed network in 
response to the driving frequency. Our work is closely 
related to the theories of stochastic resonance, ratchets 
and stochastic transport [t| but should be viewed in 
the context of synchronization. 



II. PERTURBATION THEORY FOR MARKOV 
RATE PROCESSES 

A Markov rate process over a set of states n £ 
{1, . . . , L} is described by its matrix W of transition rates 
and a master equation 

L 

Pnn (t\t ) = W nm (t)P mno {t\t ) (1) 

m— 1 



for forward-transition probabilities P nno (t\to) to be in 
state n at time t when the system is in state no at a pre- 
vious time to- The matrix of transition rates is W nm > 
for m ^ n and J2 n W n m = 0, which guarantees conser- 
vation (J2 n Pnm = 1) an d nonnegativity P nm > of the 
probability at all times. Also, because of this property 
one eigenvalue of W(t) is always zero. If this eigenvalue is 
non-degenerate, which is the case when the Markov chain 
is strongly connected, all other eigenvalues have negative 
real part and the transition probability is asymptotically 
independent of the initial conditions: 

Pnn o (*l*0 -00) = p n (t). (2) 

The perturbation theory is applicable when W(i) is 
weakly dependent on time, i.e., for W(t) = W° + eV(t) 
and sufficiently small e. From J2 n W nm (t) = for all e 
follows J2 n Vnm(t) — at all times. The time-dependent 
perturbation ansatz for the problem 

P(t\t ) = [W°+eV(t)]P(t|t ) (3) 
is to expand P(t |£o) as a series in powers of e 

P(*|*o) = P (0) (^o)+£P (1) (^o)+e 2 P (2) (^o) + ... ■ (4) 

Inserting this ansatz in Eq. §3§ and sorting by powers of 
e the dynamics of any P^ l \t\to) is determined as 

P (l \t\t ) = W°pW(t|io) + VVpV-VWto), (5) 

which, given V(t) and P ( '- 1) (t|t ), is an inhomogencous, 
linear ordinary differential equation that can be solved 
iteratively for each perturbation order. 

The mean frequency of a single stochastic oscillator is 
proportional to a sum of directed and time-averaged 
flows (J nm (t)) t = (W nm (t)p m {t)) t - In general, we write 

ft = 27r^e„ m (J„ m (t)) t , (6) 

where 0„ m £ { — 1,0,1} defines a Poincare section and 
the direction in which the Poincare section is crossed for 
each transition (FigJTJa). Note that f2 depends on the 
choice of the Poincare section, which can be arbitrary. 
Therefore, in general, is not equal to one of the 
relaxation frequencies of the system. In fact, for the two 
state Markov model of stochastic resonance [Tl|, which 
has the same form as Eq. JT]) , relaxation is not oscillatory 
at all, whereas f2 ^ 0. 

An alternative way to define mean frequency f2 
and, in addition, phase diffusion constant D, of the 
stochastic oscillator is to lift the stochastic jump process 
to a periodic lattice with W n +L,m+L = W nm . It is 
necessary to decide whether each transition is forward 
or backward (FigdJa). Then the mean frequency and the 
phase diffusion constant arc derived from the asymptotic 
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behavior of the mean and variance of n(t), for instance, 
with respect to distribution P nno (t + r\t) as 



2tt l 
il = — lim -E [n(t + r) - n ] 



2tt\ 2 1 

— lim — Var [n(t + t) 

_L / r-voo 2r 



(7) 



no 



where E and Var denote expected value and variance, 
respectively. Because of asymptotic independence from 
the initial conditions, it is possible to perform a suitable 
average over the initial conditions no and t. The 
perturbation theory for this rate process over an infinite, 
periodic lattice requires the Bloch functions of the 
infinite, periodic difference operator obtained by the lift 
of W°. We will only derive expressions for our simplest 
example, which is a one-dimensional ring with discrete 
rotational symmetry. For this example, the Bloch waves 
are plain harmonics. Dimensionless quantity c = Q/D is 
called the Peclet number and describes the coherence of 
a stochastic oscillator. Its value indicates the number of 
rotations for an ensemble of stochastic oscillators until 
coherence is lost. 

Note that, in the linear order of the perturbation, 
all time averages depend linearly on the time average 
of V(t). Therefore, without loss of generality, we can 
absorb this time average into matrix W° and assume 
(V(i)) t = henceforth. Nonlinear effects arc first 
observed in the time averages of second-order pertur- 
bation terms P^ 2 ', fi( 2 ), and of the transition 
probabilities and other quantities. We further assume 
that the perturbation is of the form 



V(i) = Ez(t)+R*z*(t), 



(8) 



where complex driving signal z(t) has rotational symme- 
try such that 



(z(t)) t = (z(t + T)z(t)) t = 0, and 
(z*(t + T)z(t)) t = e A \ 



(9) 



Equations (J8j) and © include the case of harmonic driv- 
ing z(t) = exp(— iuj-\t) where A = iu)\ and driving with 
a continuous or discrete stochastic signal with finite cor- 
relation time |Rc[A]| _1 < oo and relaxation frequency 
Im[A] = uj\. Complex matrix H assigns relative strength 
and phase to the perturbation of each transition. Thus 
we refer to H as the driving protocol. Second-order shift 
il^ of the mean frequency is a real- valued function ex- 
pressed as a sum of products between single entries of 
H* and H. A Hermitian operator corresponds to this 
quadratic form over the space of complex matrices, and 
wc define- the notation 



ft (2) (H* , H) = 2tt [v(E* , H) + v* (H* , H)] 



(10) 



where the asterisk denotes the complex conjugate. The 
perturbation expansion of the transition probabilities 



and related quantities are written in terms of the left and 
right eigenvectors v( fc ) and u^ fc ) of W°, respectively, and 
corresponding eigenvalues A&, i.e., W°u( fc ) = AfcU^'^ and 
vW+W = A^v^^, where f denotes the complex con- 
jugate transpose. We assume normalization v^ fc ^u^ fc > = 
Skk', completeness J2 k u' k 'v™t = Ilxl and conventions 
A = 0, p (0) = u (0) and v (0) = 1, where Ilxl is the iden- 
tity matrix of size L and 1 = (1,...,1) T . In Appendix 
A, we derive 



i/(A*,B)=-£e nm J>W \A 



k^O 



v (fe)tBu(°) 
nm A fc + A 



(A fc / + A) A fc 



(11) 



-w° V 

fc'#0 



for the bilinear form v(k , B) in Eq. ([10]). Equations (|T0|) 
and (|11[) capture the essence of frequency adjustment to- 
ward the frequency of the driving signal (Sec. IIV[) and 
also demonstrate the resonance nature of synchroniza- 
tion. Denominators (Afc + A) in the sums in Eq. (|11[) are 
minimized when Im[A] = — Im[Afc], i.e., when the driv- 
ing frequency matches the relaxation frequency of some 
mode k. 



III. OPTIMIZATION 

This section describes a corollary technique to de- 
termine the complex driving protocol H that optimizes 
the second-order perturbation response of the mean 
frequency. In Sec. IIVD1 we apply this technique to a 
Markov rate jump process on a random network. 

Being able to express the time-averaged second-order 
responses as a Hermitian form in terms of eigenvalues 
and eigenfunctions of the unperturbed transition matrix, 
we can formulate an optimization problem that can be 
solved using linear algebra. Consider a quadratic form 



/(x) = x^Fx 

with x e C N and a Hermitian matrix or operator F ' 
Linear constraints may be given as 



(12) 



F. 



0, 



(13) 



where the column vectors of P £ C Nxl with I < N 
can, without loss of generality, be assumed to be or- 
thonormal, i.e., P^P — l; x ;- The nullspace of P is 
spanned by another set of orthonormal vectors given by 
the column vectors of matrix Q 6 £Nx{N-i) g^jj fag^ 
Q^Q = l(Ar-z)x(AT-z) an d P^Q = 0. The purpose of op- 
timization is to maximize or minimize /(x) subjected to 
the linear constraints in Eq. (| X3[) and the constraint on 
a cost function given by a positive definite matrix S as 



x T Sx= 1. 



(14) 
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This problem is related to the variational problem in 
quantum mechanics to minimize or maximize 



$ f M$ subject to ||$|| 2 = 1 



(15) 



for some Hermitian operator M. It can be easily shown 
that the solution $ opt is the normalized eigenfunction of 
M with the largest or smallest eigenvalue. To put our 
optimization problem in this form, we need to project 
into the nullspace of P and apply a similarity transforma- 
tion that changes the constraint in Eq. (fT4|) on the cost 
function into the normalization condition on $. This is 
achieved by the following rules: 



S = Q f SQ , 

M = S- 1 / 2 Q t FQE _1/2 , 
x = QST 1/2 $ . 



(16) 



Note that with positive definite S, projected matrix E is 
also positive definite so that £ -1 / 2 is well defined and 
Hermitian. 

When applied to the problem of finding an opti- 
mal driving protocol, vector x E C L represents complex 
matrices H £ C LxL . The linear constraints can be 
as simple as requiring H nomo = for impossible 
transitions from state mo to hq or transitions that 
cannot be perturbed. More subtle constraints can be 
assigned on the phase relations. For example, setting 
H nimo = e tn Hn 2 m describes periodic switching in 
preference between two possible target states that 
are reached from state mo, which is a very simple 
model for intersections with a traffic light. Optimizing 
the frequency and phase relations between traffic lights 
in a road network is a classic problem in transport theory. 

The quadratic form /(x) = x^Fx that is to be op- 
timized can, for instance, be second-order frequency 
shift ft( 2 )(H*,H). By using the explicit expressions 

given in Eqs. (JTDJ) and (fTTj) . we can also optimize 
at a fixed driving frequency uj\. A driving protocol 
that maximizes the response to changes in the driving 
frequency can be said to have good synchronizing 
properties. 



IV. EXAMPLES 

A. Jump process on a ring 

In this section, we present the results for the mean fre- 
quency and the phase diffusion constant for a sequential 
change of states with discrete rotational or translational 
symmetry of transition rates W°. We consider L states 
that are visited in sequence with forward transition 
rates W® +1 n = L(l + 7) and backward transition rates 
W°_i n — £7 (FigJ5^,). We directly consider the lift of 



(a) 



(b)_js, 



bias wL„=L 



diffusion w d J ±t 

perturbation 

V„ +ln - A(t)cos[k„n-<p(t) 



FIG. 2. (Color online) (a) Transition rates for a biased jump 
process on a one-dimensional periodic lattice. The transition 
rate is divided into a time-independent forward bias, diffu- 
sion part, and perturbation depending on time and state, (b) 
Transition rates for a forward jump process (7 = 0) on a 
two-dimensional periodic lattice with attractive coupling in 
two directions. The speed of the transitions is indicated by 
different emphasis of the arrows. 



the finite state oscillator to the periodic lattice with 
translational symmetry W®+i m +i = W® m . The time 
scale is chosen such that the mean frequency of the 
unperturbed oscillator is loq = 2tt. The transition rates 
can be split as W° = W fl + W df into forward jump 
process W fl and unbiased diffusion part W df (see Fig(2^i). 
As discussed in Sec. HVCi these two rate processes give 
rise to the deterministic part and the diffusion of the 
stochastic oscillator in the continuum limit. 

Suppose that the perturbations are given only to 
the forward bias in the form of a traveling wave 



V„ 



W* m A(t) cos(fc ™ -</>(*)), 



(17) 



where phase <p(t) and amplitude A(t) may be determin- 
istic or stochastic. In terms of Eq. flSJ), we have 



ikom 



z(t) = A(t)e 



i<t>(t) 



(18) 



and z(t) must satisfy Eq. ©. Because of the trans- 
lational symmetry of the system the eigenfunctions of 



(fc) 



2tTU. 



(fc) 



W" and W at are harmonics given by v. 
cxp (ikn) with wave numbers k € (— 7r, 7r]. The eigenval- 
ues are given by 



W fl u (fe) = A£u« = L {er lk - 1) u 
W df u« - A d V fc ) = L 7 (e ik + e 
wV fe > = X k u k = (A« + \f )u« 



(A') 



—ik 



(19) 



Furthermore, the action of the perturbation on an eigen- 
mode creates two eigenmodes of wave numbers k ± fco, 
i.e., 



Hu( fe > 
ITuW 



2 A k+k 

1 A fl 1 

2 A fe-fc u 



,(fc+M 



(fc-fco) 



(20) 
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In Appendix B, by using the perturbation expansion of 
the characteristic function for the random process, we 
determine mean frequency f2, phase diffusion constant D 
and Peclet number c = Q/D up to the second order in 
perturbation strength e as 



UJ 



2tt, D 



2^ 



o 



(l + 2 7 ), c =u Q /D , (21) 



D 



= 1 4 



A k 



A -fco 



A A- 



A* 



(22) 



4(1 + 27) 



,^(Ag +i2 7 Im[Ag ]) 



(A.„+A)2 
2L 



A fen + A 



c.c. 



(23) 



i- = l 

co 



= 1 



0(2) D (2) 



1 



2(1 + 27) 



(l-7)A[} n +X 
A fcn + A 



Aj K+*2 7 ImK o ]) 
(A. +A)2 



+ c.c. 



(24) 



where c.c. denotes the complex conjugated terms. The 
six parameters of this model are length of the ring L 
(or alternatively the time scale of the jump process on 
the infinite lattice), diffusion parameter 7, wave num- 
ber ko of the driving traveling wave, dissipation rate 
— Re [A] and frequency Im[A] of the driving signal, and 
the perturbation strength e. If we drive deterministi- 
cally with z(t) = cxp(— iwit) exponent A = iuj\ is purely 
imaginary. The shifts of the mean frequency and the 
Peclet number for this case are shown in Figs[3Ji and 
[3b, respectively. If the stochastic oscillator is coupled 
to the mean field from a finite ensemble or a single other 
stochastic oscillator, A will have a negative real part that 
adds to the negative real part of Afc . Thus the reso- 
nance effect decreases because the absolute values of the 
denominators in Eqs. (|221) - ([24|) increase. For instance, 
with z(t) = exp(— ikini(t)), depending on another un- 
perturbed oscillator with stochastic jump process n\{t), 
the exponent of the autocorrelation function is simply 
A = \-k t [see Eq. (|92|) in Appendix B]. Driving with 
another stochastic oscillator of identical length, i.e., with 
similar stochasticity but variable mean frequency, the fre- 
quency shift is weaker, and we observe no enhancement 
in the Peclet number (see Figs. [3J; and|3ji). 



B. Mutually coupled oscillators 

Unlike an externally driven stochastic oscillator, 
system n(t) = (no(i), ni(i)) of two mutually coupled 
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FIG. 3. (Color online) Relative frequency change (left col- 
umn) and change in the Peclet number (right column) of a 
stochastic biased jump process on a ring in (a) and (b) for 
the case of harmonic driving of the forward transition rates 
with A = iuj\ , in (c) and (d) for driving with another stochas- 
tic oscillator of the same length L and A* = A df + A fl ^i/27r, 
and in (e) and (f) for mutually coupled oscillators of same 
length L, k° = (2ir/L,-2w/L) and /3 = tt/2. Numerical 
simulations were performed using a generalized Gillespie al- 
gorithm for time-dependent transition rates. Averages were 
taken over 1000 runs for 1000 time units each with the pa- 
rameters 7 = 0.1, e = 0.5. The inset in subfigure (a) applies 
to all subfigures. 



stochastic oscillators is autonomous. As we pointed 
out, it is permitted to take an average over the initial 
conditions in Eqs. (JT)). So far, in the time-dependent per- 
turbation theory, averaging over time has considerably 
simplified the equations. Now, the perturbation does 
not depend on time, i.e., z(t) = 1 in Eq. ([S]). Instead, 
we average f > n o +n,n o (' r |0) over the initial conditions 
no (see Appendix C). We have verified our results 
against Monte Carlo simulations (FigsOJ; and [3J) of 
the coupled rate processes performed with the usual 
Gillespie algorithm for time independent transition rates. 

We again consider the lift of the combined system 
n = (no, n\) = rioe + nie 1 to an infinite periodic lattice 
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(see Fig|2jD). The unperturbed transition rates may be 



given by isotropic diffusion part W df and forward bias 
W fl with different values in the directions e° and e 1 : 

w° = w fl + w df , 



(25) 



Coupling shall be given by matrix 

Vn+e°,n = L COs(k°n + 0) , 



V n+e \ n = ^Lcos(k°n-^). 



(26) 



Term k°n denotes the inner product of n with wave vec- 
tor k°. In contrast to external driving, here phase shift 
j3 does have physical effects. For k° = (27r£ — 1 , —2itL~ 1 ) 
and /? = 7r/2 coupling is attractive, whereas for f3 = 
cosine coupling inhibits synchronization. The former 
case is presented in Figsj3k and[3f. 



turbed system are harmonics = {^) 2 u^ = e lkn . 



Due to symmetry, the cigenfunctions of the unper 

turbed sy; 
We define 

\df _ \df | \df 
A k — A k + A fci' 



A k (P) = h + ^ X ki e > 

A k = A«(0) + A df , 



(27) 



where and A df are given by Eq. p^|) and k = (k , ki). 
In Appendix C we determine the frequency, phase diffu- 
sion constant and Peclet number of the first oscillator to 
the second-order in perturbation strength e as 



n 

CJ 



= 1 - £ z 



1 



A u o 



-k° 



D 
D 



1 



4(1 + 2 7 ) 



*k»H3) (A2o+i2 7 Im 



(28) 



A 1.0 



X 2 



C 
CO 



2(A« +L) + A«o(-/3) 



A k o 



A« -7A5o(-/3)+i 



(29) 



2(1 + 27) 



A k o 



\ 2 



(30) 



In contrast to the driven system, frequency wi of the 
second oscillator enters both the denominator and 



numerator of the expressions via A k o and A k0 (— P), 
respectively. As a result, the second-order perturbation 
terms do not vanish in the limit lj\ — > oo, since the 
phases of the fast and the slow oscillator remain corre- 
lated. Another remarkable difference is the possibility 
of an increase in coherence for two mutually coupled, 
identical or weakly nonidentical stochastic oscillators 

(FigEF). 

While the time-dependent perturbation theory can 
be used to approximate the frequency and the phase 
diffusion constant under driving by another stochastic 
oscillator, the transition rates in the combined system 
(no(t), ni(t)) are not time dependent, irrespective of 
mutual or unidirectional coupling. The mean frequency 
in the combined autonomous system can, in principle, 
be obtained nonperturbativcly from the stationary prob- 
ability distribution, i.e., the explicit zero eigenvector 
of the combined time-independent matrix of transition 
rates. In [l2j this was performed for a master-slave 
pair of two-state oscillators, and frequency locking 
was observed for strong coupling. The phase diffusion 
constant cannot be obtained in such a simple way. In 
fl3| a method has been recently found to determine 
the phase diffusion constant in a periodically driven 
system from the cyclostationary solution of an extended, 
periodically driven, linear ordinary differential equation. 
In that case, the lift to the periodic lattice and the 
generating function approach is not necessary. It would 
be worth investigating whether this method can be 
generalized to a larger number of states, nonsequential 
transition networks, or stochastic driving. 



C. Continuum limit 

The purpose of this section is to illustrate that our 
theory is accurate in the limit of weak coupling and 
strong noise and complements the Kramers rate theory 
in the other asymptotic regime of strong coupling and 
weak noise. 

In the limit L — > oo, the biased jump process on a 
ring network with harmonic driving corresponds to a 
continuous, weakly perturbed limit-cycle oscillator in 
the Kuramoto phase approximation 



i9 = 27r(l+ecos(i9-o;it)) + y/2Do^(t) (31) 

with delta correlated white noise (£(£)£(£')) = S(t — t'). 
To obtain a finite phase diffusion constant, j/L must be 
kept constant such that, in the limit L — > oo, 



47T 2 7 

L 



D . 



(32) 



The second-order perturbation terms for the frequency 
shift, the phase diffusion constant and the Peclet number 
in the continuum limit arc given in Appendix B. However, 
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FIG. 4. (Color online) Plot of tilted potential (bold, blue line) 
in the stochastic Adler equation (|36|l . The local minimum at 
<po also determines the two adjacent local maxima at — ipo 
and 2ir — ipo. The potential difference over the left barrier is 
— 2U(ipo) and 2irA — 2U(ipo) over the right barrier. 



exact nonpcrturbativc expressions for these quantities 
can be found. Substituting tp = -d—uiit and A = [oj\— 2tt) 
in Eq. ([ST]) , we obtain the stochastic Adler equation 



tp = -A + 2tt£ cos ip + 



(33) 



for which the mean frequency and the phase diffusion 
constant are known in terms of special integrals fl4l . [i~5| . 
Here, we will use the continuum limit to compare the 
different regimes in which our perturbation theory and 
the Kramers theory are valid. Noting that we use a time 
scale such that loq = 2n, we introduce dimensionless pa- 
rameters x = A/ Do and y = 2ire/Do- The continuum 
limit of the second-order perturbation frequency shift in 
Eq. (|2"21 is given in Appendix B [Eq. (|96b[) ]. In terms of 
x and y, Eq. (|96b[) can be rewritten as 



LU 



D 



y_ 

2 



1 



(34) 



The derivative of this quantity with respect to x is max- 
imal at x = with 



9. 



uj 



Do 



2 



(35) 



A value of this derivative larger than unity is unphysi- 
cal, because then the frequency response is stronger than 
the change in the driving frequency. Therefore, we say 
y = \/2 is an upper bound for the region in which our 
perturbation theory is valid. On the other hand, the 
Kramers approximation gives the rate of phase slips in 
the Arnold tongue region for small noise strength Do or 
strong coupling. Equation (|33|) can be written as a gra- 
dient system 



-U'{<p) + V2D~ !;(t), 
(A • tp — 27resiniy9) . 



(36) 



The tilted sine potential is shown in Fig(4] It has a lo- 
cal minimum at ip = arccos(x/y) > and two adja- 
cent maxima at tp\ = —ipo and tp2 = 2-k — ipo- With 
U min = U(ipo) and U max = U(tpi) or U max = U(tp 2 ), the 
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FIG. 5. (Color online) Rescaled frequency shift (Q—uo)/ Do of 
the periodically driven drift- and diffusion process on a circle 

(a) as a function of x — A/Do and y = 2ne/Do. Values x < 
and x > correspond to negative (blue) and positive (red) 
frequency shifts, respectively (color bar aligned on top). The 
dashed lines mark the border of the Arnold tongue region of 
phase locking in the limit Do — s> 0. The Kramers rate theory 
is valid between the yellow (light, solid) curves that indicate 
the condition ([/ max — C^min) /Do — 1, whereas below the blue 
(dark, solid) line for y < \/2, our perturbation theory will 
give physically sensible results for all frequency differences, 
in particular, at resonance point A = 0. Above the blue 
(dark, solid) line, the perturbation theory may be still give 
good results sufficiently far away from the resonance point. 

(b) Frequency shift as a function of y with x = 1.0 fixed 
(solid curve). The two dashed curves indicate the second- 
order perturbation approximation for low values of y and the 
Kramers rate approximation for high y. The exact frequency 
shift was determined by numerical evaluation of the integrals 
in 0. 



Kramers theory requires (£/ max — Umin) /Do 3> 1- We 
can parameterize the condition ([/ max — £7 m i n )/-Do = 1 
by < ipo 5: 7r /2 and obtain 



x = ± ycosipo- (37) 



2 (sin ipo — ipo cos ipo) ' 
Kramers rate tk over a potential barrier is given (l6j by 

^ sJu;UU;Ue {U -- U — )/Do ■ (38) 



From = (6) t = ui + (ip) t and by subtracting the back- 
ward jump rates from the forward jump rates, we obtain 



Q — ujq 
Do 



x -yam(<p ) e 2 ^°- vs ' m ^ (l - e - 2lTX ) 



(39) 

Figure [S] compares the regions in the (x, y) parame- 
ter space where the Kramers approximation holds with 
that where our second-order perturbation theory is valid 
at the resonance point. Small noise or large coupling 
strengths will result in strongly nonlinear behavior near 
the resonant frequency and render the perturbation the- 
ory invalid. In contrast, the regime of the Kramers ap- 
proximation can always be reached by increasing the cou- 
pling strength, or if \x/y\ < 1, by reducing the noise 
strength. Figure [SJs shows the rescaled frequency shift 
at x = A/ Dq = 1 as a function of y = 2-nejDo- At low 
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FIG. 6. (Color online) (a) Unidirectional ring network with 
L = 30 states and E sc = 15 random shortcuts, (b) Rearrang- 
ing the nodes on the unit circle according to the argument of 
the complex eigenvector to eigenvalue Xmax = — 0.52 — 0.77i of 
maximal ratio r ma x = 1-48 between imaginary and real part 
strong circular flows in the network can be detected. For- 
ward transitions are shown in blue (solid arrows) and back- 
ward transitions in red (dashed arrows), (c) Convex hull of 
the complex eigenvalues excluding Ao = for single realiza- 
tions of directed ring networks of size L = 1000 at various 
shortcut densities a = E sc /L. Here E sc is the number of 
additional random transition channels (including duplicates). 
The dashed lines originating at zero are tangents of slope 
\r\ = 1.09 to the convex hull for a = 1.0. The dots indi- 
cate nonzero complex eigenvalues for a network realization at 
a = 2.0. The spectral gap increases with a and the slope 
of the tangent decreases, (d) Double logarithmic plot of the 
average maximum value of r = Im[A]/Re[A] in the spectrum 
versus shortcut density a. Each data point and standard devi- 
ation was determined from a sample of 1000 random networks 
of size L = 500. The dashed line in (c) is r — 1/y/a. 



coupling strengths y < the second-order perturba- 
tion theory approximates the quadratic behavior of the 
frequency shift, whereas at high coupling strengths, the 
Kramers rate theory describes the exponential deviation 
from frequency locking. Therefore, it is appropriate to 
consider the adaptation of frequency in our theory as a 
weak form of synchronization. 



D. Ring with random shortcuts 

While the jump process on a ring is a discrete 
approximation for a limit-cycle oscillator with finite 
phase diffusion, this one-dimensional geometry may not 
reflect the topology of a more complicated oscillatory or 
transport process. A way to incorporate a more complex 



FIG. 7. (Color online) Optimization of the driving protocol 
for the network shown in Fig(6b. Assigning the same rate 
and perturbation cost to each transition in (a), we plot the 
eigenvalues of the derivative of the frequency shift operator 
(thin lines) with respect to frequency u}\ of deterministic har- 
monic driving. In (b), the eigenvalues of the second-order 
frequency-shift operator (thin lines) are shown as functions 
of the driving frequency. The red (bold, light gray) and blue 
(bold, dark gray) lines in (a) and (b) show the frequency shift 
and its derivative when the driving protocol is chosen to have 
the optimal frequency response at either of the two apparent 
resonant frequencies seen in (a). 



topology is to use a network of states that are not 
visited in sequence. Here we study the effect of adding 
E sc random shortcuts to a unidirectional ring of size L 
(Figj6ji). The resulting Markov chain is strongly con- 
nected ensuring that eigenvalue Ao = is non-degenerate. 

The strength of the resonance depends on coher- 
ence ratio rfc = Im[Afc]/Re[Afc] of resonant eigenmode 
k. The corresponding eigenvector u.( fc ) and its complex 
conjugate quantify collective oscillations and are asso- 
ciated with nonzero stationary probability fluxes in the 
system. Figure [6^ shows a typical network for L = 30, 
E sc = 15, i.e., shortcut density a = 0.5, where each node 
n is located on a circle at phase d n = 2im/L. The same 

network is shown in FigO) with i? n = axg(un), where 
u( fc ) is the right eigenvector of W° for the eigenvalue 
with the maximal ratio r^. Forward transitions are 
highlighted blue (solid) and backward transitions red 
(thin, dashed). With respect to this arrangement, the 
flow is strongly biased in the forward direction. 

We investigated the dependence of the expected 
maximum value r max of the coherence ratio on shortcut 
density a = E sc / L. The eigenvalues of the transition 
matrices are localized in a region of the complex plane 
with a negative real part and are symmetric with respect 
to the line of real numbers (Fig[6j;). The slope of a line 
originating at zero and tangent to that region scales 
with some power of a in both limits a — > and a — > 00 
(FigJBJi). The decrease in r max with increasing a is 
the result of stronger mixing because of the random 
shortcuts. In the region of the topological crossover 
at a = 1 we observe that the coherence ratio is 
approximately equal to unity (see FigslHt and[BJi). 
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When the transition rates are perturbed harmoni- 
cally with frequency u>\ and driving protocol H that 
encodes the relative amplitudes and phases between 
the perturbations, the mean frequency will change by 
a value given approximately by Eqs. (fTUj) and (fTTj) in 
SecHU Using the technique described in Sec lIIII we 
optimized the driving protocol for the example shown 
in FigJBjD such that, for a fixed driving frequency <Ji, 
responsiveness of the frequency shift is maximized. 
We assumed that the perturbation cost is the same for 
every transition. We chose ui\ to be one of the apparent 
resonant frequencies in Figj7^i. Synchronization, i.e., 
the adaptation of the mean frequency to the driving 
frequency, both positive and negative, is observed in 
FigEb when these driving protocols are kept constant 
and uj\ is changed. 



V. DISCUSSION 

We have presented second-order perturbation analysis 
for the frequency of periodically or stochastically driven, 
time-continuous Markov chain models. In the case of 
a biased jump process on a ring, we also determined 
the second-order perturbation terms for the phase 
diffusion constant. This simplest discrete-state model 
of a stochastic oscillator can adapt its mean frequency 
toward the frequency of a deterministic or stochastic 
driving signal. The mechanism for this weak form 
of synchronization is a parametric resonance of the 
driving signal with an oscillatory relaxation mode of 
the stochastic process. We also showed that, when 
two identical stochastic oscillators arc attractively 
coupled, they can increase their Pcclct number, a 
measure for the coherence of oscillations, above the 
level of the uncoupled system. Furthermore, the explicit 
expressions for the frequency shift could be used to op- 
timize the perturbations in a complex transition network. 

Phase response can be viewed as a perturbation 
theory for the mean frequency of an oscillator. The tech- 
niques presented in this paper could therefore be useful 
in developing effective phase models for noise-induced or 
strongly noise-perturbed oscillations (l5j. The explicit 
perturbation expressions for the frequency shift of a 
periodically driven Markov rate process may be used to 
estimate the unperturbed transition rates. In the case of 
enzymatic reactions, it may thus be possible to examine 
the reaction kinetics by measuring the response of the 
turnover rate to changes in the frequency and wave form 
of the driving signal. 



the transition rates in the asymmetric exclusion process, 
which models transport through a cell membrane, was 
studied perturbatively. The resonance of the mean flow 
with the driving frequency in the second order in the 
perturbation strength was found, which agrees with our 
general results. 



The authors thank Lutz Schimansky-Geicr and Jun 
Ohkubo for helpful comments. 



APPENDIX A 

In this section, we will derive general expressions for 
the time-dependent transition probabilities of a driven 
Markov rate process up to the second order e 2 in pertur- 
bation strength. Transition probabilities P nno (t, to) obey 
the following equation 



ith 



P(t\t ) = [W° + sV(t)] P(t|i ) 



V(i) = Rz(t)+R*z*(t). 



(40) 



(41) 



We assume that complex driving signal z(t) is stationary 
and symmetric such that 



(z(t + T)z(t)) t = (z(t)) t = 0, and 
c(r) = (z*(t + r)z(t)) t = e AT forr>0. 



(42) 



The perturbation expansion of P(t |£q) will be given in 
terms of eigenvectors and eigenvalues of unperturbed 
transition rates W . We denote the left and right eigen- 
vectors as v( fc ) and vS k ^ , respectively, and the correspond- 
ing eigenvalues as The completeness of and orthogo- 
nality relations for the eigenvectors are ^2 k u^v^t = 1 
and v^^u^ •* = Skk', respectively. In the main text we 
also set v ( n ] = 1 for 1 < n < L. If W° is an infinite peri- 
odic difference operator, the eigenvalues are distributed 
continuously in a finite number of bands. Then, the com- 
pleteness relation is expressed as a sum over different 
bands and for each band an integral over wave numbers 
(—it, tt\. The orthogonality relation is given by the prod- 
uct of a Kronecker delta for the band number and a Dirac 
delta function for the wave numbers within the bands. 
The nondegencrate zero eigenvalue corresponding to the 
unique stationary solution of the unperturbed problem is 
denoted as Aq = 0. By inserting the perturbation ansatz 



P(t\t ) = P^(t\t ) + Y,e l P {l) (t\t ) 



(43) 



Of course, the theory presented in this paper is 
also applicable to other nonequilibrium systems that are 
reasonably well described by stochastic cycles through 
the states of a Markov chain as well as stochastic trans- 
port or growth processes. In [l8j . periodic driving of 



in Eq. (|40|) and sorting by powers of e, the dynamics of 
each perturbation term is given by a linear, inhomoge- 
neous, ordinary differential equation 

P (l) (t\to)=W P < - l Ht\to)+V(t)P {l - 1) (t\t ). (44) 
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In terms of Fourier modes 

42'(*l*o) = v«tp«(^ )u( fc '), 

q+ k , =v( fe )t H u( fc '), (45) 
q ~ kk , = v^W' 

and using the completeness relation, Eq. (|44p becomes 

42< - AfcTrW + E feW*) + 4'^. ( 46 ) 

k" 

Note that because v^H = v^H* = 0, coefficients g± , 
vanish for k = 0. We will explicitly exempt these modes 
from the sum or integral in Eq. (|46p . From the initial 
condition P nm (t\t) = P^l(t\t) = 5 nm follows 42' = 
Sio^kk' and therefore 



(47) 



All higher order perturbation terms are given by convo- 
lutions of the inhomogeneous part with an exponential 
kernel. Defining 



(48) 



the general solution to Eq. (|4*5]) for I > is 



k" 



r a-i) 



r (i-i) 
'fe"fc' 



(r) 
(r). 



(49) 



The first- and second-order perturbation terms are thus, 
respectively, 



4 1 Mt + r\t) = e^ q + k , J C t kk ,[z](r) 
+ e XkT q kk ,£iAz*](T), 



(50) 



and 



(z*(t + T) Z (t + t')) t e -( A *-V)f 'dt' 



c(r-t')e" (Afc " V) *'^' 
c(t / )e (A *'- A * )(r -*'W 



= e (V-^£ fc , fe [c](T) 



(52) 



and 



£L" [ z *£fc»fc'Mj ( T ) 

= r e (v-v)t'£^„[ c ](t') e -(A,-A fc „)f^/ (53) 

Jo 

z jCkk' jCk'k"l c ] ( r )- 



For c(t) = exp(Ar), the incomplete Laplace transforma- 
tions are solved explicitly. For short hand notation, we 
define 



9kk' = *v - A fe + A. 



Then, 



and 



(54) 



(55) 



T Cw Civic] ( t ) 



1 



9k' k" 



e W'* _ 1 



jCkk' 

x kT ( e g' k \,,T _ 1 g^,,, 



(56) 



1 



gbk' 



9t" 



k"^0 



For fc' = k, this equation takes the form 

X * T ( e 9 kk" T - 1 



*$(t + r|t) = e^ J2 ltk"%"k'C\k" UCl'Az*]} W eAfcT ilfefc (r 



SfcV 



fffefc" 



+ E ikk"ii"k' elk" [**&»*[*]] ( r ) 

fc"^0 



+ 0(z 2 



(51) 



(57) 

For fc 7^ and k" ^ 0, the dependence on the initial 
conditions is removed by taking limit r — > oo in Eqs. (|55[) 
and (l56l: 



Here 0(^ 2 ) denotes the products of the form z(t+t')z(t+ 
t") and z*(t + t')z*(t + t") which will vanish when aver- 
aged over t. The remaining nonvanishing time averages 
are expressed by incomplete Laplace transformations of 
complex autocorrelation function c(r) = (z*(t + r)z{t)) t 
or its complex conjugate. We observe 



lip e Vr £° fc a[c](r) = -^o4- 

1 



lime A -£° fefe ' (r)=*. 



(58) 



->fe'0 



A, 



The limits exist because the real parts of Afc and Afe" are 
strictly negative and the real part of A is assumed to be 
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less than or equal to zero. Terms irj^i with k' ^ vanish 
in the limit r — > oo. The explicit time dependence of the 
transition probabilities is necessary to obtain the gener- 
ating function for random process n(t). From the gen- 
erating function, we can calculate both mean frequency 
and the asymptotic phase diffusion constant, as shown in 
Appendix B for the biased jump process on a ring. In 
the general case, the mean frequency can also be calcu- 
lated from the asymptotic mean probability flows in the 
system. Let Q„m G { — 1,0,1} define for each transition 
whether it crosses a Poincare section and in which di- 
rection. Mean frequency oj is proportional to the time 
averaged flow over this Poincare section times 27r, 



and 



ft = 2-7T ^ Q nm (Jnm(t)) t 



+ e 2 (wZ m p$(t) + V nm (t)p£\t) 



(59) 



(60) 



The time averages in the linear order of e are zero be- 
cause of Eq. (|42[) . We arc now able to express the 
mean frequency in the second order of e 2 by using 
Eqs. (|5H 1) .(15T1 ) .([55 ]) ,(|56 ]) . and (J5SJ. We have 



V nm (t)p$(t) 

-V/) (h 

fe#0 



r (1) 



H 



run 



r (1) 

T k0 



,(0)* 



fe#0 \ 9 k 9 0k / 



(61) 



as well as 



kk'^0 



( k )W° ( lkk'lk'0 , ^kk'Ik'O 
Afcffofc' ^kgQ k , 



(0), 




,(0)* 



(62) 



Shift fi' 2 -' in the mean frequency to the order e 2 as a 
function of the complex driving protocol H is given in 
terms of the bilinear form z/(A*,B) 



(63) 



ft (2) (H* , H) = 2tt [i/(H* , H) + i/* (H* , H)] 



The bilinear form i/(A ,B) is obtained from Eqs. ([54] 
and ([55 ]) -([52 ) as 



v (fc)tB u (°) 
1 A fe +A 



v (fc)t A *u( fc ')v( fc ')tBu(°) , 

(^Ta^ I 



If exponent A depends in some way on frequency uj\ of 
complex driving signal z(t), the responsiveness of the fre- 
quency shift to oj\ is calculated from 



d ul v(A*,B) = d Ul Aj2®n m J2 U ™ V 0) A n 



nm k^0 

-w° V 

r r nm / j 



v (fc)tB u (°' 
1 (A, + A) 2 



I • (65) 



k'^0 



(Afc' + A)' A fc y 



APPENDIX B 



We consider a biased jump process on a ring of states 

L(l + 7) and 
Z/y. The stationary 
probability flow in this system is J° + i „ = n — 

W° n+1 )/L = 1, i.e., the time scale is chosen such that 
the mean period is one and the mean frequency is ujq = 
2n. We split the transition rates into an unbiased jump 

Lj 

The perturbation shall be given 
a s a traveling wave to the forward transitions 



with forward transition rates W® +1 n 
backward transition rates W° n+1 



process with forward and backward rates W£±i n 
and bias W^ +l n = L 



V nm {t) = W* m A(t) cos(fc m - (f>{t)). 



(66) 



Amplitude A(t) and phase 4>(t) of the perturbation spec- 
ify complex driving signal z(t) = A(t) exp(— i(j>(t)), which 
is supposed to be symmetric and posses a discrete or con- 
tinuous rotational symmetry. Then, the eigenvectors of 
the unperturbed system are simple harmonics and the 
perturbation couples a mode k to its neighboring modes 
k ± fco . The mean frequency in the second order in per- 
turbation strength can be calculated from the expres- 
sions derived in Appendix A. Here we will use a different 
approach and derive the frequency shift as well as the 
asymptotic phase diffusion constant from the generating 
function of the biased jump process on an infinite lattice. 
In this case, eigenmodes 



,(fc) 



27™i fc) 



Akn 



k G ( — 71, 7r] 



(67) 



form a complete and orthogonal set with 

v {k) *u {k) dk = S nm and 



(68) 



v«t u ( fe '> = ]T «W*u^ =5(k-k'). 



The eigenvalues follow from 



V w° 



u {k) 

nm ^m 



Xk.ui k) 



= L(l + 1 ){u^_ l -u^)+L 1 (u^ +l ~u^ 



[£(1 + 7) 



11, 
I) 



(69) 



, i k 



,0) 
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They can be divided into diffusion and flow parts as 
Afe = A? + Xf, 



A fe > ^k 

Xl = L (e-^ 1) , 
Xf = L-{ (e lk + er lk - 2) 



(70) 2tt 



Perturbation V = Hz(i) + H*z*(<) applied to mode u( fe ) 
creates two neighboring modes 

u (fc±fc ) as 



H*n( fc ) - ix a ii( fe - fe o) 



(71) 



so that 



g +,=vWtHu( fe ') = iA^(fc-fc'-fco), 



v (fe)t H * u (fe') = I A fl 5(fe _ fc ' + fco) . 



(72) 



The mean frequency and the phase diffusion constant are 
defined as 

2tt 1 
n = — lim —E\n(t + r) - n(t)], 

/2tt\ 2 1 (?3) 

D = f -T I lim 7^Var([n(< + r) - 

\ L J T^roc 2t 

Taking the long time limit makes to and D independent 
from the initial conditions. Therefore, we can start with 
a localized distribution at n(t) = and average over time, 
i.e., study at the asymptotic behavior of the moments of 
the time-averaged transition probability (P n o(t + r\t)) t . 
These moments arc conveniently calculated from the 
characteristic function 

<P T {x) = E [e ixn ] = £ e ixn (P n0 (t + r\t)) t (74) 



as 

E[n] = -i4>'(0), and Var[n] = -0"(O) + /2 (O). (75) 

The perturbation expansion of the time-averaged transi- 
tion probabilities is 

(Pno(t + r\t)) t = P$(r) + e 2 (p$(i + r|i)) t ■ (76) 

In Appendix A, we have already derived the Fourier 
modes of the transition probabilities up to the second 
order of e 2 . 



(77) 



E 



e-<((^ t+ e 2 (^)jvr*dkdki 



Here, we have used 

1 oo 

£ e i( fe+ x)n = ^ = ^(fc.fc/). (78) 



n— — oo 



From Eq. (|5"Tj) . we see that the second-order perturba- 
tion term of average Fourier coefficient (tt^^ ^ contains 

products of the form <ikk"1k"k'i an< ^ with Eq. (|T2")) . it fol- 
lows that these terms are nonzero only if k = k' and 
k" = k =p kg . All integrals over the Fourier modes are 
thus evaluated explicitly, and we find 



= 4 6 



■= + fc r — 1 



A I A 

j — x, — x-\-ko \ x,— x+ko 



(79) 



i. e A-xr A -x-k A -x I e 



o A * I o A * 

^— x,— x — ko \ ^—x,—x—ko 



Inserting Eq. (|T9")) into Eq. (|77|) . we find that the gener- 
ating function has the form 



x) = e A -» r (1 + e 2 h(x)) 



(80) 



To emphasize the structure of h(x) and take the deriva- 
tive, we introduce some more shorthand notations 



h+(x) 



1 A« 



x+k 



4 9 A . t : 

•3 — X,— X-\-KQ 



h+{x) = A« 



K(x) 



9 ~x,—x-\-ko 
J- A -a:-fco 

4 <? A * 

^ — x, — x — ko 



(81) 



h (x) = X 



—x A* 

^— a: — x — fco 



so that 



h(x) =hf (x)fi2 {x) + h l (x)h 2 (x) 

-(h+(x)+h-(x))X ti _ x T. 



(82) 



The necessary derivatives of the eigenvalues are 



dx ~ x+k ° 



a 



dx - x+k <> 



dx 



-A df 

2 *-x+k 



-2 7 Im[A«J, 



(83) 



x=0 
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and with Eqs. fl80|) and ([75]). it follows that 
E[n] = Lt- e 2 iti(0) and 

Var[n] = L(l + 2 7 )r - e 2 h"(0). 



(84) 



Functions hf(x) and all their derivatives at zero vanish 
when divided by r in the limit r — > oo 

1 <f 



lim 

r^oo r ax 



0, (7i = 0,l,...). (85) 



Only derivatives of hf (x)\^_ x t at zero remain in that 
limit. With Eq. (j8l"j) and recalling that g£ k , = Afc/ — Afc + 
A, we have 



hf(0) 



1 A fcn 



4 A fcn +A' 



(86) 



hno) = i-, 



*t+L A« o (A« o +z2 7 Im[Afl o ]) 



Afc n + A 



(Afc +A)2 



and hj (0) = ftf(0)* and ftf'(O) = /if'(0)*. Finally, we 
can collect all the terms necessary to write explicit ex- 
pressions for the mean phase velocity, the phase diffusion 
constant, and their ratio the Peclet number 



2tt 1 2ir 1 

tt = — lim -£[n] = 2tt - e 2 — lim -ih'(Q) 

= 2 . + £ 2 ^(^(0) + / l r(0))z^A fi ; 



(87) 



2tt 1 - e A 



1 



A fco 



A 11 



A feo + A A_fc + A 



and 



'2tt\ 1 
D = | — lim — Varfnl 
L ) r^oo 2r L J 



1 



so that D Q = (1 + 2 7 )27r 2 /L and 



D , 1 

= 1 + £ 

£>o L(l + 2 7 ) 



«f (0) S A», 



2=0 



+ c.c. 



1 



£C=0 

'+rn\ - h+i 



(1 + 2 7 ) 
1 

4(1 + 2 7 ) 



2iti+(0)~h+{0) + c.c. 

Ag (Ag o +z2 7 Im[Ag Q ]) 
(A fco +A)2 



3Ag +2L 
A fco + A 



+ c.c. 



(89) 



The relative Peclet number is 

"ft( 2 > £>( 2 >1 



c 

CO 



1 



w 



= 1 



1 



2(1 + 2 7 ) 



(1-7)A[ 1 + L 



A 



feo 



A 



(90) 



X t ( A fc + l2 Tlm [Ag ]) 



+ c.c. 



(Afc„+A)2 

where c.c. denotes the complex conjugated terms. 

From the generating function Eq. ([74]) also follows the 
autocorrelation of the driving signal if it is taken from 
another stochastic oscillator as z(t) = exp(— ifcni(t)) : 

(z*(t + r\t)z{t)) t = J2 e ikn (Pno(t + r\t)) t , (91) 



(92) 



i.e., A = A_fc = XI in Eq. (|42|) . Because the time scale 
of the driving signal can be different from the perturbed 
stochastic oscillator, eigenvalue A_fc is not necessarily 
given by Eq. ([7DJ). 

The continuum limit can be taken if L, 7 — > oo 
and 



and with Eqs. (jjjj and <J80J), we find 

(z*(t + T\t)z(t)) t = m=e X - kT , 



2ir 2 

Do = — (l + 2 7 ) 



(93) 



is kept constant. For harmonic driving, angle $ = 2ixnjL 
in the continuum limit evolves according to 



i? = 2tt(1 + wit)) + V5i^£(f)i (94) 

where £(t) is white noise with (£(t)£(t')) = S(t-t'). After 
substituting 93 = $ — Wit, we obtain the stochastic Adlcr 
equation 



99 = 2tt - wi + 2tt£ cos(<^) + ^/2D £,{t) (95) 

for which the mean freq uency and the phase diffusion 
constant are known [3 EHl- With ko = 2n/L and 
A = uji — 2ir, the eigenvalues and the perturbation terms 
become 





■> - 


i2?r, 


n 


-> 1 


+ £ 2 7T 


UJQ 






D 


■» 1 




AT 


+ £ 2 



A 



df 



£ 2 + A 2 ' 
(2n) 2 (A 2 - Dq) 
(D 2 + A 2 ) 2 ^ Z?q + A 2 



2^ 2 



(96a) 
(96b) 

(96c) 
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APPENDIX C 



and 



Here we consider two stochastic jump processes no(t) 
and ni(i) on a two-dimensional, periodic lattice. While 
we have used the time-dependent perturbation theory in 
the case of external driving, the coupled system n = 
(^Oj^i) = n o e ° + nie 1 is autonomous. Nevertheless, 
because the time-independent perturbation theory is a 
special case of the time dependent perturbation theory, 
we can set z(t) = 1 and follow the derivations in Ap- 
pendix B with a slight modification. The time averages 
have to be replaced with a suitable average over initial 
state n° = n(0). For this, we introduce unitary trans- 
lation operator T n o with T^oPTnO = P n o +n n o +m . 

J nm 

Transition probabilities P(r) = P(i + r\t) are now inde- 
pendent of t, but Fourier modes 

7r kk '(n ,r) = v( k )t T | i0 P(r)T n ou( k ') 

and 

^(n°)=vWtTt HT n ou( k '), 

? kk '(n°)=v (k)t T| 1 oH*T n ou( k ') 

depend explicitly on the initial conditions. Due to 
translational symmetry, transition matrix W° and the 
translation operator T n o commute. 



1 



(97) 
(98) 



?kW kn ^A k (/3)£(k-k'-k°), 
9 kk '=^ k ° n ° V^A k (-/^(k-k' + k°). 



(103) 



We see that the average over the initial conditions in 
Eqs. (f5U|) and (f5T|) removes all terms linear in q^ k , and 

quadratic terms ? kk «<7 k » k ' anc ^ 9kk"9k"k'- The generating 
function is written as 



(x)=^ e - n <P n o +n , n o(T)) n0 



(104) 



and from Eqs. d5TJ) , d5TJ) , and (|103[) and setting A = 0, we 
have 



f (2) 



dk' 



(105) 



l c A_ xr Afl x+k0 (-/3)A fl x (/3) / e (^ +k o-A- x )r_ 1 

4 A_ x+k o — A_x V A_ x+k o — A_ x 

1 x T A fl x _ k o(/3)A fl x (-/3) ( e &— 



r 



A 



-x- k ° 



A_, 



A 



-x- k ° 



A_, 



Again, we divide the transition rates into isotropic 
diffusion part W d and a bias in the forward directions. 



w° = w df + w fl , 

df 



^„± e o,n - W? ±eKn = L 7 , 



(99) 



The oscillators are coupled weakly through their forward 
jump rates with the perturbation 

Vn+e°,n = L COs(k°n + /3) , 

w lr „o n, (100) 



£cos(k°n-/3). 

Z7T 



Here k°n denotes the inner product. The eigen- 
functions of the unperturbed system are harmonics 



(k) _ 

11 6 


lkn . With k = 


(ko,fa 




= L (e- lk - 1) 


5 


Af 


= L (e lk + e- 1 


*-2) 


A k 


_ xdf , x df 
- *k„ + A fc l5 




A k (/3) 


_ \& , Ul \A 

- A k„ + ^fc, 




A k 


= A k (0) + A k f . 





(101) 



Then 



W°u( k ) = A k u( k ), 
Hu( k ) = i e ^A k+k0 (/3)u( k + k °), 

H*u( k ' = i e -^A k _ k o(-/3)u( k - k °) 



(102) 



As in the previous section, the generating function is of 
the form 



(x) 



1 + £ 2 /l(x)) 



(106) 



and the mean and variance of the position of the first 
oscillator are given, respectively, by 



E[n ] = -i <9x <?K x )lx=o ' 
Var[n ] - - ^(x)| x=0 + (^ o 0(x)| x=o ) s 

With the definitions 

1 A fl x+k0 (-/3) 



4(x) = A fl x (/3) 

KM 



4 A_ x+k o - A_ x 



A- x + k o — A_> 
1 A« x _ k0 (/3) 



ft 2 -(x) = A« x (-/3) 



4 A_ x _ k o — A_ x ' 



A- x - k o — A_ x 
function /i(x) can be written as 

ft(x) = /i+(x)/i+(x) + }q{x)}q (x) 

-(^(x)A fl x (/?)+Mx)A fl x R3)) 



(107) 



(108) 



(109) 



Again, only the derivatives of terms /ij = (x)A^ x (±/3) re- 
main in the limit r — > oo. The derivatives of the eigen- 
values with respect to the first component of x are the 
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same as those for the single oscillator 

9x A^. x+k o(/3)| x=0 = faT^-xo+k* 



f) X df I 

o x A -x+k°| 



x df 



=° dx X0+k " 



Xn=0 



X =0 



(110) 



given by Eq. ([55)1 in Appendix B. Thus, functions /if (x) 
and their derivatives at zero are 



h+(0) 



1 A* (-/3) 



4 A k o 
9 Xo h+ (x)| x=Q 



(111) 



A«„+L A k o(-/3)(A«„+i2 7 I 



in 



X 11 



A 



k° 



A^.0 



so that 



D , 2 

7T = 1 + e 

^0 



2a co /if(x)a ro A_ x (/3) 



= l + e 



L(l + 2 7 ) 

+ / l +(x)a 2 A_ x (/3)+c.c. 
1 



x=0 



(1 + 2 7 ) 
1 



2id Xo h^(x) - /if (x) + c.c 

AjoC-g) (Ago+z2 7 Im[Ago]) 



x=0 



X 2 



4(1 + 27) 

2(A£ +L) + A k oR3) 



A k o 



+ c.c. 



(114) 



as well as h~(0) = h+(0)* and d Xo h-(0) = d Xo h+(0)*. 
The mean frequency and phase diffusion constant of the 
first oscillator are then determined as 

2tt 1 2n 1 

fi = ^— lim — E[ra ] = 27r — e 2 — lim — id xo h'(x)\ a 

L t— >-oo J" X/ r^oo T 

= 2tt + e 2 ^ [/if (x)3 X0 A« x (/3) + /h (x)a K0 A fl x (-/3)] | x=Q 



and for the Peclet number 



= l + S 2 



( n 


D 






l 





2(1 + 27) 



'A^o-7A k „(-/3)+i 



A k o 



2tt 1 - e 1 



A k o(-/3) , A« k0 (/3) 



A k o 



\_ k 



A ko(-/?) (AfeO+i2 7 Ii 



x 11 



(112) 



X 2 



+ c.c. 



(115) 



£> = — lim — Var[n] 

T— >00 2T 



'2tt 

27T 2 

L 



(113) 



(1 + 27) 
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